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A SIMPLE AND EFFECTIVE HIGH-ORDER SHOCK-CAPTURING 
LIMITER FOR DISCONTINUOUS GALERKIN METHODS 

SCOTT A. MOE*, JAMES A. ROSSMANITHt, AND DAVID C. SEAL* 


Abstract. The discontinuous Galerkin (DG) finite element method when applied to hyperbolic 
conservation laws requires the use of shock-capturing limiters in order to suppress unphysical oscil¬ 
lations near large solution gradients. In this work we develop a novel shock-capturing limiter that 
combines key ideas from the limiter of Barth and Jespersen [AIAA-89-0366 (1989)] and the maxi¬ 
mum principle preserving (MPP) framework of Zhang and Shu [Proc. R. Soc. A, 467 (2011), pp. 
2752-2776]. The limiting strategy is based on traversing the mesh element-by-element in order to 
(1) find local upper and lower bounds on user-defined variables by sampling these variables on neigh¬ 
boring elements, and (2) to then enforce these local bounds by minimally damping the high-order 
corrections. The main advantages of this limiting strategy is that it is simple to implement, effective 
at shock capturing, and retains high-order accuracy of the solution in smooth regimes. The resulting 
numerical scheme is applied to several standard numerical tests in both one and two-dimensions and 
on both Cartesian and unstructured grids. These tests are used as benchmarks to verify and assess 
the accuracy and robustness of the method. 


1. Introduction. 


1.1. Some background. Godunov’s theorem states that a Zmear numerical 
method when applied to the one-dimensional, linear, advection equation: 


( 1 . 1 ) 


dq 

di 


^=0 

dx 


q{t = 0,x) = qo{x), 


where a; S M, t € K>o, and q{t, x) : K>o x M —>• K, is monotone (i.e., does not generate 
new extrema) only if it is at most first-order accurate. In practice, schemes that 
are not monotone will exhibit large spurious oscillations near discontinuities in the 
solution. The analog of this result for more complicated hyperbolic equations (i.e., 
nonlinear systems in multiple dimensions) is that linear numerical methods with a 
formal accuracy higher than first-order will create spurious oscillations that can lead 
to catastrophic numerical instabilities (e.g., negative densities and/or pressures). As 
in the linear scalar case, these instabilities are especially pronounced if the solution 
becomes discontinuous (i.e., shock waves and contact discontinuities). The additional 
challenge with nonlinear hyperbolic equations is that generic initial conditions tend 
to produce shocks in finite time. 

The principal remedy for these numerical instabilities in high-order methods is to 
introduce a shock-capturing limiter, such limiters render an otherwise linear method 
nonlinear, thereby circumventing the restrictions of Godunov’s theorem. Typically, 
limiters are designed to execute two tasks: (1) to monitor the solution quality (at least 
in some simplified sense); and, when necessary, (2) to reduce the influence of high- 
order corrections. The development of limiters for hyperbolic equations was pioneered 
by Bram van Leer in a series of papers starting with [40{[^ . Additional fundamental 
contributions came from Harten [9 10 , Sweby 


38 , and Tadmor 39 who developed 


mathematical precise limiting strategies based on the principle that the numerical 
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scheme should be total variation non-increasing. A good review of these limiters can 
be found in Chapter 6 of LeVeque [22| . In more recent work, many contributions have 
been made to generalize these limiting strategies to high-order methods, including 


to weighted non-oscillatory (WENO) schemes (e.g., see Liu et al. 23 and Jiang and 


Shu ) and discontinuous Galerkin schemes (e.g., for a review see Qiu and Shu 
and Dumbser et al. (!)• 

In this work we focus on limiters for discontinuous Galerkin (DG) methods 


31 


Limiting strategies for DG methods can generally be broken down into two key 
steps: 

Step 1. Detect troubled-cells (i.e., cells that contain large gradients); 

Step 2. Apply a limiter to each troubled cell that is detected (i.e., in some way reduce 
the local gradients or higher derivative analogs thereof without changing the 
average value in that cell). 

In some limiting strategies these two steps are executed with completely different 
techniques (e.g., see [^), while in other approaches the detection and limiting happen 
in all in one fell swoop (e.g., the moment limiter of Krivodonova 16 ). In either case, 
the end goal of most limiters is to design a method that simultaneously accomplishes 
two often conflicting objectives: (1) reduce unphysical oscillations in the presence of 
shocks and (2) retain high-order accuracy in smooth regimes. 

The techniques used for DG methods in the limiting step (Step 2) can be subdi¬ 
vided into three broad classes of limiters, which we briefly describe below. 

1.1.1. Moment or slope methods. These limiters have their roots in the 
development of classical second-order finite volume methods. The key ingredients are 
minimum modulus evaluations that compare the slope (or higher derivative analogs) 
in the current cell to values computed from neighboring cells. Such methods are still 
widely used (e.g., see [l2p3p6p8H20p4] ). There are two well-known difficulties with 
these limiters: (1) they tend to degrade the of order of accuracy in smooth extrema; 
and (2) they are difficult to generalize to unstructured meshes. 

1.1.2. Hybrid weighted essentially non-oscillatory methods. Methods in 
this category are motivated by the success of the weighted essentially non-oscillatory 


(WENO) method in shock-capturing. The classical schemes 27,29,30,48-51 rede¬ 
fine the polynomial representation of the solution inside each element by considering 
neighboring elements and matching cell averages. The WENO reconstruction, which is 
automatically mass conservative, is then typically applied in order to redefine polyno¬ 
mial representations inside each element. More recently, the so-called Hermite WENO 


schemes have been used in an effort to reduce the size of the computational stencil 24 


A chief criticism of these methods has been the expense of their application; although 
recent efforts exist to reduce their computational cost 51 . 


1.1.3. Artificial viscosity methods. 


idea of adding a small amount of artificial viscosity to hyperbolic equations 42 


This method is based on the classical 

The 


and further developed in 


modern incarnation of this limiting strategy was pioneered by Persson and Peraire 26 

EEl 


These methods typically work by using a smoothness 


indicator to first detect troubled cells, and then directly discretizing a diffusive term 
that has been artificially added to the equations in an effort to smooth-out oscillations. 
One drawback of using these limiters with explicit time-stepping schemes is that due 
to the diffusion operator, the maximum stable time-step is At = 0{h/M^) (i.e., 
the viscosity parameter is typically 0{h)), where h is the grid spacing and Md is 
the highest polynomial order of the basis polynomials, rather than the typical At = 
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0{h/Miy) for hyperbolic systems. 


Other variations such as minimum entropy satisfying limiters 47 have also begun 


to be investigated. A complete description of what options exist for limiting DG 
solutions is beyond the scope of this work; we direct the interested reader to for a 
recent and extensive review of the literature. 

1.2. Brief summary of current work. The purpose of this work is to de¬ 
velop a novel limiter for discontinuous Galerkin methods that has all of the following 
properties: 

• Provable high-order accuracy in smooth regimes; 

• Robust reduction in spurious oscillations in the presence of shocks; 

• Computationally efficient (i.e., only uses neighboring cells and no character¬ 
istic decompositions); 

• Easy to implement in existing DG codes; 

• Automatic extension to multidimensional settings including Cartesian and 
unstructured grids; and 

• Automatic incorporation of positivity preserving limiters. 

The limiter presented in this work can be viewed as a novel extension of two separate 
methods. First, it can be thought of as extending the finite volume Barth-Jespersen 
limiter to the discontinuous Galerkin framework. Second, it can be viewed as an 
extension of the modern maximum principle preserving DG schemes developed by 
Zhang and Shu 46 . In particular, the limiter is based on finding for each element 


local upper and lower bounds using only nearest neighbor values, and then limiting 
the higher-order modes in such a way that that mass conservation is automatic and 
such that the limited solution remains within the predicted local bounds. Roughly 
speaking, this limiter can be viewed as belonging to the moment or slope methods, 
although it has important differences with the current methods in this class, including 
no need to for characteristic decompositions, no need for hierarchical highest-to-lowest 
moment limiting, and easy extensibility to unstructured meshes. 

The remainder of this paper is organized as follows. In Section[^we briefly review 
the DG method. The basic limiting strategy for a one-dimensional scalar equation is 
outlined in Section]^ In Section|^we describe how the proposed limiter obtains high- 
order accuracy. Extension to multiple dimensions is shown in Section and to general 
systems of hyperbolic conservation laws in Section where we use the compressible 
Euler equations as the prototypical example. In Sectionj^we describe how positivity¬ 
preserving ideas can be incorporated. The numerical method is validated on several 
standard cases for the compressible Euler equations in Section Gonclusions are 
made in Section [9l 

1.3. Implementation details. All of the numerical examples carried out in 
this work are implemented in the software package DoGPack 


32 . This software 


package is open-source and all the code used in this work can be freely downloaded 
from the web. 

2. DG-FEM framework. In order to define the proposed limiter, we first pro¬ 
vide a brief description of the discontinuous Galerkin (DG) method, which serves to 
define the notation used throughout the remainder of this work. We refer the inter¬ 
ested reader to Idlll 111361 and references therein for further details of the DG method. 


Gonsider a conservation law of the form 


( 2 . 1 ) 


9.t 


V-F(g)=0, 

3 


in n C 







with appropriate boundary conditions, where q{t, x) : K>o X !->• are the 

conserved variables, F(q) : i—)• jg function, d is the spatial 

dimension, and Me is the number of equations in the system. We assume that the 
system is hyperbolic, meaning that the flux Jacobian matrix, n • F.^, is diagonalizable 
with real eigenvalues for all q in the domain of interest and for all ||n|| = 1. 

Let C be a polygonal domain with boundary 917. The domain 17 is dis¬ 
cretized via a finite set of non-overlapping elements, %, such that 17 = Let 

pMn denote the set of polynomials from to M with maximal polynomial 

degree Md. Let denote the broken finite element space on the mesh: 

(2.2) W'* := G : zc'*|.^^G , \J% G T^} , 

where h is the grid spacing. The above expression means that G has Me 
components, each of which when restricted to some element 71 is a polynomial of 
degree at most Md and no continuity is assumed across element edges (or faces in 
3D). 

The approximate solution on each element 7) at time 7 = i" is of the form 


(2.3) 


<7^r,x(0) 


Ti 


Ml(Md) 

i: <3‘'”v''(«). 

e=i 


where Ml is the number of Legendre polynomials and (^) : i—)■ M are the 

Legendre polynomials defined on the reference element To in terms of the reference 
coordinates ^ G Tq. The Legendre polynomials are orthonormal with respect to the 
following inner product: 


(2.4) 


1 






'To 


1 if k = £, 
0 if k^t 


We note that independent of h, d, Md, and the type of element, the lowest order 
Legendre polynomial is always = 1. This makes the first Legendre coefficient the 
cell average: 


(2.5) =\hIro 


=: C- 


To obtain the semi-discrete DG method, equation (2.1) is multiplied by the test 
function the resulting equation is integrated over an element 7), and integrations- 
by-part are performed on the flux term. The result is the following system of ODEs: 


( 2 . 6 ) 


m\. 


'Ti 




where the numerical fluxes T must be determined via (approximate) Riemann solvers. 
In this work we make use of the Rusanov 33 numerical flux, which is often referred 


to as the local Lax-Friedrichs (LLF) flux in the literature. 


Finally, in order to obtain a Runge-Kutta DG method, equation (2.6) is discretized 


in time via some Runge-Kutta method. In this work we make use of strong stability¬ 
preserving (SSP) Runge-Kutta methods, such as those found in Gottlieb, Shu, and 
Tadmor and Ketcheson 


15 
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3. Proposed limiter: The scalar case. The basic procedure we follow to limit 
the solution on a single element Ti is given by the following sequence of steps: 

Step 0. Define a set of points Xi that will be used to approximate cell maximum and 
minimum values. In the current work, we use Gaussian quadrature nodes and 
augment them with corner points and quadrature points along the element 
boundaries (i.e., edge Gaussian quadrature points). 

Step 1. For each mesh element, 7^, we compute an approximate maximum and min¬ 
imum: 


qM, ■= max q^{x) 
xex 


Ti 


and qm '■= min < q {x) 
xGx 


Ti 


Step 2. We consider the set N-j-i of all neighbors of Ti, excluding % itself, and com¬ 
pute an approximate upper and lower bound: 


(3.1) 

(3.2) 


Mi := max qt -|- a{h), max { qu, } ? , 
rrii := min qi - a{h), min {qm, } > • 


The scalar function a(h) > 0 is a tolerance function that will be described 
shortly. The most aggressive limiter we consider sets a = 0. 

Step 3. We define 


(3.3) 


OMi '■ = 


- qi 

QMi - Qi 


and 9rri, := 


rrii - qi 
Qrrii Qi 


where 0 < 4>{y) < 1 is a cutoff function. In this work, we use the function 
(3.4) </)(?/):= min I^,l| , 


which is motivated by the work of 25 . The form of this function becomes 
important in the multidimensional setting, which will be described shortly. 
Step 4. We define the rescaling parameter as 

(3.5) 0* := min{1, 6»m., 6 »mJ ■ 

Step 5. Finally, we rescale the approximate solution on the element 77 as 


(3.6) 


T{x) qi + 0^ ( q’^{x) 


Ti 




Before delving into the finer details of the method, a few remarks are in order. 
Remark 1. The choice of quadrature points in Step 1 is not unique. Indeed, 
the purpose of this step is to construct an approximate upper and lower bound for the 
solution. We originally tested this limiter by using cell averages, but found that those 
solutions tended to be quite diffusive near shocks. Directly sampling neighboring cells 
(that are potentially oscillating) allows us to retain sharper features in the solution. 

Remark 2. The presence of a non-zero value of a in Step 3 is required to obtain 
high-order accuracy at local extrema. This issue will be elaborated upon in Section IH 
Remark 3. When the proposed scheme is compared to the recent maximum 
principle preserving (MPP) limiters significant differences come into play in 

Step 2 and Step 3. 
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• For the MPP methods, Step 2 is replaced by fixed a priori global bounds, 
whereas we follow the Barth-Jespersen idea and use this step to estimate 
local upper and lower bounds for the solution. 

• In Step 3, the MPP limiters use hard upper and lower bounds for the solution 
that are known a priori, whereas we use the result from Step 2 to estimate 
these. 

• In the MPP literature, authors normally set 4'{y) = min{l, 2 /}, whereas we 
find that the form of this function is important when pushing to multiple 
dimensions. In Section!^ we elaborate on this idea. 

Additionally, the newly proposed limiter is designed to capture shocks, whereas the 
original MPP limiter for DG methods was solely designed to preserve the maximum 
principle property of hyperbolic problems. In their work, additional limiting was 
necessary in order to suppress unphysical oscillations. 

Remark 4. The rescaling in Step 6 does not effect the total mass. 

The consequence of this is that the limiter will be automatically mass conservative, 
which is an important property for hyperbolic solvers. This can be observed by simply 
integrating equation (3.6) over a single element. 


4. On the importance of a for retaining high-order accuracy. In this 
section we describe the importance of a for obtaining a method that is genuinely high- 
order accurate. Many limiters exhibit the so-called “clipping” phenomenon, where 
smooth extrema are cut off because the limiter turns on and damps out these values. 
Recent efforts have been put forth to mitigate this clipping phenomenon for other 
discretizations, including the piecewise parabolic method (PPM) [^. 

In this section, we will argue that high-order accuracy will be retained provide that 
the scalar function a (/i) : K>o —>■ R>o vanishes slower than Oih?). This observation 
was inspired by recent work on extensions of the Barth-Jespersen limiter to finite 
volume methods [25| , and we demonstrate the performance of our limiter by testing 
it on a smooth solution in ID. 

In one dimension we can break down any smooth function on a bounded domain 
into monotonic regions and finitely many points of local extrema. We begin with a 
discussion of monotonic regions, and then turn to isolated local extrema. 


4.1. Performance in monotonic regions. In regions where the exact solution 
is monotonic, the proposed limiter does not effect the asymptotic convergence rate 
independent of the choice of a > 0. Without loss of generality, we consider a region 
with the solution is non-decreasing. After after enough refinement, the exact solution 
satisfies qi-i < Pi < Qi+i where qi refers to the solution value on the z*^-grid cell in a 
monotonic region. The algorithm then correctly selects qi+i as an upper bound, qi-i 
as a lower bound, and therefore does not overreach by over limiting the solution on 
the grid cell. A formal proof of the order of accuracy would follow that already 
presented in 46 , but would start with the fact that our algorithm over-estimates the 
upper bound and under-estimates the lower bound for the exact solution. 


4.2. Performance in regions near smooth extrema. In this subsection, we 
describe the importance of using a non-zero value of a to maintain genuine high-order 
accuracy. If a is chosen to be too large, then oscillations in a given cell will never 
be clamped down, whereas if a vanishes too quickly, then smooth extrema will be 
clipped. Our goal is to find necessary conditions for maintaining high-order accuracy 
at smooth extrema. 
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To begin, we consider a Taylor expansion of a smooth function g : K —>■ K at an 
extrema x = ^o- 

q{x) = q{^o) + q'{^o){x - Co) + ^q”iCo)ix - + ■ • • ■ 

Given that q'{Co) = 0, we observe that asymptotically |g(x) — (;(^o)| = 0{h?) if 
|x —^o| = 0{h). This leads us to consider functions of the form a{h) = 0{h^), 
where r < 2 in order to maintain 


\q{x)-qiCo)\^0{h'^) <\a{h)\ for all \x-Co\=0{h). 

Note that in the 2D Cartesian case: h = max (Ax, Ay), where Ax and Ay are the 
mesh widths in each coordinate direction. 

That is, as long as a{h) goes to zero slow enough, then our limiter will eventually 
turn off and not clip smooth extrema. However, if we were to set a = 0, then our 
limiter will induce clipping at smooth extrema, and yet if a does not vanish as h —>■ 0, 
then our method would permit problematic oscillations to persist. 

The choice of a is not unique, and if desired, could be tuned for any given problem. 
The trade-off is that if a is too large, then true oscillations will only be damped 
provided the mesh is fine enough. Our choice of a{h) = is consistent with 

the recommended switches that turn off similar Barth-Jespersen type limiters for 
finite volume schemes 25 . Extensive numerical tests indicate this choice finds a good 


balance between the need to limit near shocks, while being soft enough to allow high- 
order accuracy on coarse grids. In the case of smooth problems, our limited 
solution will always converge to the unlimited solution even near local 
extrema. We illustrate this point with a simple case study. 

4.3. On the importance of a: A simple ID case study. We analyze the 


performance of the proposed limiter on the linear advection equation (1.1) on x G [0,1] 
with periodic boundary conditions and with the initial condition 


(4.1) 


y(x) = 


j cos® 

(x—0.5)71 \ 

. 016 J 

lo 



if |x-0.5| < 0.08, 
otherwise. 


In Table [2 we present a convergence study for various grid resolutions after advecting 
the solution through one full period (e.g., solve to t = 1). For this example we use the 
10-stage SSP fourth-order Runge-Kutta scheme (SSPRK4) developed by Ketcheson 
15 and hold a constant CFL number of 0.4. We compare the results for the unlimited 
solution, the solution with a = 0, and the solution with a{h) = 50h^'^. 

We observe that after enough refinement, the limited and unlimited solutions 
agree with each other, whereas for very coarse grids, the order of accuracy is slightly 
degraded. It is only after the limiter “turns off” that we we dial in on high-order 
accuracy. This cutoff point is a function of the curvature of the exact solution at the 
local extrema. This example serves to illustrate that after enough refinement, there 
is no difference between the limited and unlimited solutions. 


5. On the importance of (p for multidimensional cases. An important 
distinction between the single and multidimensional cases is that in higher dimensions, 
shocks are no longer isolated to single cells. To illustrate this, consider the following 
initial condition: 


(5.1) 


q{x,y) 


1, if X < 0, 

0, otherwise. 
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Mesh 

No Limiter 

Order 

a = 0 

Order 

a = 50/ii-3 

Order 

5 

3.78 X 10-^^ 

— 

8.73 X lO-'J^ 

— 

3.78 X 10-oi 

— 

8 

5.61 X 10-°i 

-0.84 

9.00 X 10-01 

-0.07 

5.61 X 10-01 

-0.84 

14 

2.32 X lO"*^! 

1.58 

8.22 X 10-01 

0.16 

2.34 X 10-01 

1.56 

24 

5.38 X 10-°2 

2.71 

7.04 X 10-01 

0.29 

1.74 X 10-01 

0.56 

42 

4.24 X 10"^^ 

4.54 

4.78 X 10-01 

0.69 

7.82 X lO-o^ 

1.43 

73 

3.65 X 10"°^ 

4.44 

2.24 X 10-01 

1.38 

2.19 X 10-02 

2.30 

127 

3.89 X 10"'^'^ 

4.04 

8.84 X lO-o^ 

1.68 

2.73 X 10-03 

3.76 

222 

4.10 X 10-°6 

4.03 

2.97 X 10-02 

1.95 

4.10 X 10-06 

11.64 

388 

4.38 X 10"^^ 

4.01 

9.37 X 10-00 

2.07 

4.38 X 10-07 

4.01 

679 

4.65 X 10"°® 

4.01 

2.76 X 10-03 

2.19 

4.65 X lO-o® 

4.01 


Table 1 


Convergence for the ID advection problem in Section \4-3\ All errors are L 2 norm errors. We 
see that a nonzero value of a is required to allow the limiting to completely turn ojf when the solution 
is properly resolved. 


Now, consider a Cartesian mesh with cells of width Ax and height Ay: 


(5.2) Q, := 




1 

*+2 


Ax 


j - ^ j Ay, (j + ^ 


Ay 


i,j e Z. 


After projecting the initial conditions onto with Md > 1 (i.e., see definition 


(2.2)), every cell Coj Vj will exhibit Gibbs phenomena; this is because the initial 
discontinuity bisects the cells Coj Vj. We now consider the effect of applying our 
limiter to this solution. We focus on what happens to the solution on cell Coo. 

In Step 2 of the proposed limiter, we compute upper and lower bounds of the 
solution using nearest neighbors, which in this case satisfy 


max 




That is, for these initial conditions, the predicted upper bound for cell Coo will not 
clamp down on any oscillations present in this cell, unless we are careful about deciding 
on the correct structure of (j). 

The important observation to make here is that it is the ratio of the deviation from 


cell averages (the value 9 in equation (3.3) of Step 3) that is the relevant quantity to 
study. For this problem, it is identically equal to 1 for every mesh element. Indeed, 
what this tells us is that it is necessary for (/>(!) < 1 if we will have any hope of 
clamping down on the unavoidable oscillations. We observe that even with ^(1) < 1, 
the oscillations are not completely suppressed after a single time step; however, after 
several time-steps the artificial oscillations decay exponentially to zero. In practice, 
we observe that this is sufficient to construct robust simulations. 

As a simple demonstration of the importance of this choice of (/>, we consider 
linear advection in two-dimensions on a Cartesian grid. Consider the scalar advection 
equation 

dq dq dq 
dt~^ dx~^ dy 

with double-periodic boundary conditions and discontinuous initial conditions: 


(5.3) 


+ r^ + ^=0, for (x,y) G [0,1] X [0,1], 


(5.4) 


qo{x,y) = 


1 if a: G [0.3, 0.7], 
0 otherwise. 





























q(t x,y)att=1 



q{t, X, 0.5) att = 1 












< 




0 0.25 0.5 0.75 1 

X 


(b) ID slice of 0(y) = 



Fig. 1. Square wave test case for scalar advection with double periodic boundary conditions. 
The solution in Panels (a) and (b) is computed using 0(y) = min{l,y}, while the solution in Panels 
(c) and (d) is computed using (i){y) = min{l, 3 ^}. The simulation is run to a final time oft = 1, at 
which point the solution should return to the initial conditions. The choice in (a) and (b) satisfies 
0(1) = 1, whereas the choice in (c) and (d) satisfies 0(1) < 1. The first choice permits oscillations to 
remain in the solution, whereas the second choice eventually suppresses the unphysical oscillations. 
The remainder of the simulations used in this work will use the function 0(y) = min{l, 3 ^}. 


In Figure we show results for a grid of size 60 x 60 for a fourth-order (P 3 ) method, 
comparing the choices of (p{y) = min{l, y} and 4>{y) = min{l, j^}- This result shows 
that the second choice is necessary for pushing our results to multiple dimensions. 
Other choices exist, but one property we find important is that </)(!) < 1. In the 
remainder of our numerical simulations, we choose 4’{y) = min{l, in all the 2D 
examples. 


6. Extensions to systems of equations. As a prototypical example of a hy¬ 
perbolic system of equations, we consider the compressible Euler equations: 


( 6 . 1 ) 


p 


pu 

pu 

-b V- 

puu + pi 

_£ _ 


u {£ + p)_ 


^ = ^Pl|u|P + 


P 

7-1’ 


where 7 is the gas constant. The conserved and primitive variables are 


( 6 . 2 ) q := (^p, pu^, pu^, pu ^and w := (^p,u^,u^,u^,p) , 

respectively. In these expressions p is the mass density, pu := {pu^, pu^, pu^) is the 
momentum density, £ is the total energy density, u := (u^, u^) is the fluid velocity, 

and p is the pressure. 
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The simplest implementation of the proposed limiter for this system would be to 
apply the scalar limiting procedure to each conserved variable independently. Indeed, 
our initial tests employed this approach. However, after extensive testing for this 
system, we discovered that using the primitive variables to define one high-order 
damping parameter 9 yields better results. An example is provided in Section [6.1| to 
verify this claim. 

The extension of the proposed limiter to systems of equations can be summarized 
as follows: 

Step 0. Select a set of variables to use in the bounds checking: w. In fluid dynamics, 
primitive variables are the preferred choice due to their Galilean invariance. 
Select a set of points Xi that will be used to approximate cell maximum and 
minimum values. In this work we always select: corners, internal, and edge 
Gaussian quadrature points. 

Step 1. For each element i and each component (. oi w compute: 


(6.3) 


wX[. := max < w^{q (x)) 


Ti 


and . 


min 

xexi 





Step 2. For each element i and each component i oi w compute an approximate 
upper and lower bound over the set of neighbors, Nj-. (excluding the current 
cell?)): 


(6.4) 

(6.5) 


max < 

wf + aih), max i Wm \ 


{ jeNr^ 1 J 

min < 

Wi - a{h), min 1 \ 

1 

jeNr^ 1 ^ J J 


where tcf is the cell average of over cell 7). Nj-.. In the Gartesian grid case, 
we define Nj-. to be the set of cells that share a common edge with 7). This 
provides good results in general, and has the nice property that it does not 
extend the effective stencil size. For unstructured grids, we define Nj-. as the 
set of all cells that share a common node with 7). If this choice is replaced by 
elements that share common edges only, then the results are more diffusive. 

Step 3. For each element i compute: 


( 6 . 6 ) 9m := min 



Step 4. For each element i compute: 


(6.7) 6 'j:=min{l, 9mi, ^'mJ- 

Step 5. For each element i and each component £ oi q compute 


( 6 . 8 ) 






Note that we are using the values of the variable w to determine the amount of limiting 
that needs to be done (i.e., size of 9i), but we are actually applying the limiter to the 
conserved variables, thereby retaining mass conservation at the discrete level. 
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6.1. Primitive vs. conservative variables for computing 6: A ID ex¬ 
ample. In order to test the limiter performance using either primitive or conserva¬ 
tive variables to determine 9i, we consider the one-dimensional shock-entropy prob¬ 
lem found in 23 35,37, which highlights the interplay between shock-capturing and 
feature-preservation. The initial conditions for this problem are given by 


{p, u^,p) 


(3.857143,2.629369,10.3333) 
(1 -|- esin(5a;), 0,1) 


X < —4, 
X > —4. 


We follow the common practice of setting e = 0.2 and discretize the domain [—5, 5] 
into a total of 200 cells. The final time for this simulation is t = 1.8. This test prob¬ 
lem involves a shockwave interacting with a smooth entropy wave. As this interaction 
occurs, the result is the formation of highly oscillatory waves after the shock, where 
high-order accuracy should produce a benefit, yet limiters often smear out these os¬ 
cillations. Ideally a limiter should be able to pick out and limit only non-physical 
oscillations, while simultaneously capturing the oscillatory post shock features. 

The results of two sets of simulations, one where 9i is computed from sampling 
conserved variables and one from sampling primitive variables, is shown in Figure 
For each version of the limiter, we highlight the effect of the choice of a that permits 
subcell oscillations to persist. Additionally, we show the full results of this problem 
with the primitive variables used to compute 9i and a = 500/i^-® in Figure These 
results clearly show the advantage of using primitive variables in defining 9i and give 
a range of a values that seem to give good results. 

In every simulation past this point we use primitive variables to define the rescal¬ 
ing parameter and we take a = 500h^ ®, which seems to balance the need to suppress 
unphysical oscillations and still retain high-order accuracy in smooth regions. The 
only exception to this is the example in Section [8.3[ where for most of the computa¬ 
tional domain we do use the primitive variables to calculate 9i, but there is a small 
region around the corner of the backward-facing step where we switch to a more 
aggressive limiter based on conservative variables. 

6.2. On the importance of a: A simple 2D case study. The convergence 
properties shown in Section [T^ should still be observable in 2D and for systems. To 
verify this we consider a 2D Euler example that has a smooth exact solution. This 
example involves an ideal gas with 7 = 1.4, on (x,i/) € [0,2] x [0,2] with double- 
periodic boundary conditions and the following initial conditions: 


(6.9) 


(^p,u^,u^,p) = (1 -I- 0 . 2 sin( 7 r(a: -I- ?/)), 0.7, 0.3, l.O). 


In the exact solution u^, u^, and p remain constant for all space and time, while the 
density, p, is advected with the constant fluid velocity given by {ui,U 2 ) = (0.7,0.3). 
The numerical L 2 errors after one revolution (i.e., time t = 2) are shown in Tab le [2 
To produce these results we used the same time-stepping scheme as in Section 4.3 
We again observe that we must have a nonzero value of a to asymptotically match 
the unlimited solution’s convergence even for this very simple test problem. 

7. Positivity preservation. Given the framework of the limiter, the inclusion 
of the recent and popular positivity-preserving methods developed by Zhang and 
Shu |46| is essentially automatic. For example, to guarantee positivity of a scalar 
function that has a positive mean, we simply replace the computation of the rescaling 


parameter in equation (3.5) of Step 4 of the algorithm with an additional term to 
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^"Conservative Variables 
O Primitive Variables 
• Reference 

Fig. 2. Shock-entropy problem with 4 points plotted per cell. In these panels we highlight the 
benefit of using the primitive variables to conduct the limiting at negligible additional cost. To high¬ 
light the study, we compare different results that are constructing by modifying the subcell tolerance 
parameter ex. We observe that selecting a relatively large value of a combined with primitive vari¬ 
ables is our best choice, given that the conserved variables tend to introduce additional oscillations 
post shock in smooth regimes. 


take the minimum over: 


I5 ^mi 5 ^Mi 5 


Qi 
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6i — min 


Qi Qmi 















5 


4 

3 

2 

1 

-5 - 3.75 - 2.5 - 1.25 0 1.25 2.5 3.75 5 

Fig. 3. The full plot of the shock-entropy wave problem with a = and primitive variable 

limiting. This a value seems to give a good compromise of resolution and limiting. We will use 
a = 500/i1'5 in our future calculations unless otherwise noted. 


Density with a=1000.0 

- 

Ax'"" 



.^Reference 
• DG-Primitive Variables 


Mesh 

No Limiter 

Order 

a = 0 

Order 

a = 0.5h^® 

Order 

10 

1.81 X 10-°5 

— 

4.24 X 10^°2 

— 

1.07 X 10”°2 

— 

15 

3.39 X 

4.97 

3.12 X 10"^^^ 

0.92 

6.07 X lO"*^® 

1.68 

22 

6.87 X 10"°^ 

4.48 

2.23 X 10-“2 

0.94 

2.12 X 10"°® 

2.95 

33 

1.33 X 10"^^ 

4.05 

1.75 X 10"“^ 

0.59 

8.56 X lO"*^® 

7.91 

50 

2.50 X 10"°® 

4.36 

1.12 X 10-°2 

1.17 

2.50 X 10"°® 

21.25 

75 

4.91 X IQ-'^y 

4.02 

7.84 X lO-'J® 

0.88 

4.91 X lO-'^y 

4.02 

113 

9.52 X 10"^'^ 

4.15 

5.11 X 10"“® 

1.08 

9.52 X 10"^'^ 

4.15 


Table 2 


Convergence for the 2D Euler equation problem from Section \6.2\ All errors are L 2 norm errors. 
We see that with a nonzero value of a, we eventually match the unlimited problem’s convergence 
rate, but that without a nonzero value of a, the convergence is very poor. In fact when the other 
methods are showing their theoretical asymptotic convergence rates the method using a = 0 is still 
showing first order convergence (it should be able to achieve second order convergence). 


Extensions to the compressible Euler equations to maintain positive pressure require 
expanding the unknown variables as a quadratic function of 6. The details of this can 
be found in the work of Zhang and Shu 


46 


8. Additional numerical results. In this section we apply our numerical 
scheme to several additional standard numerical test cases. We include one and two- 
dimensional test cases on both Cartesian and unstructured grids. These benchmark 
test cases are used to verify the accuracy and robustness of the proposed method. 


8.1. 2D Riemann problems. First, we consider two multidimensional Rie- 
mann problem test cases. These problems have been used extensively as benchmark 
test cases for other limiters [6 17 


The first example, RPl, is a well known example that can be found in 34 


and [^. We divide the domain [—0.5, 0.5] x [—0.5,0.5] into four rectangular quadrants 
all of which meet at the point (0.3, 0.3). We number the quadrants starting in the 
upper left hand corner, and continue in a counterclockwise manner. For example, the 
upper left hand corner is Quadrant 1, and the upper right hand corner is Quadrant 
4. The initial conditions for RPl are defined in Tabled 

The second two-dimensional Riemann problem we consider is called RP3 [6 17 


The domain [—0.5, 0.5] x [—0.5,0.5] is divided into four equal area quadrants, and the 
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Quadrant 

P 



P 

P 



P 

1 

0.5323 

1.2060 

0.0000 

0.3000 

2.00 

0.75 

0.50 

1.00 

2 

0.1380 

1.2060 

1.2060 

0.0290 

1.00 

-0.75 

0.50 

1.00 

3 

0.5323 

0.0000 

1.2060 

0.3000 

3.00 

-0.75 

-0.50 

1.00 

4 

1.5000 

0.0000 

0.0000 

1.5000 

1.00 

0.75 

-0.50 

0.00 


Table 3 


Initial conditions for RPl (columns 2-5) and RP3 (columns 6-9). 


numbering of the quadrants is also defined in a counterclockwise fashion starting at 
the upper left hand corner. The initial conditions for RP3 are defined in Table 

We run both of these examples on a grid with 400 x 400 elements and use third- 
order P 2 elements. Time-stepping was done with the classical third-order SSPRK3 
method using a CFL number of 0.1. Both examples are run using a Rusanov (local 
Lax Friedrichs (LLF)) flux [^. Results are plotted in Figure]^ which are schlieren 
plots that show the magnitude of the gradient of the solution density |Vp|. We show 
two figures, one with 1 point per cell and the other with three points per direction 
per cell. We make this comparison because the limiter permits tiny oscillations to 
remain in the solution, since it ignores all variations smaller than a certain threshold 
(controlled by the parameter a). As the mesh is refined, a also shrinks, and therefore 
the oscillations will vanish with mesh refinement. We note that these oscillations are 
not visible in the plots showing one point per cell. 


8.2. Double Mach reflection. For our next two-dimensional example, we run 
the classical double Mach reflection test problem (e.g., see Qiu and Shu 30 ). In this 
test problem a Mach 10 shock approaches a wedge with an angle of 30 degrees. In 
order to simplify the implementation of this problem on a Cartesian grid, it is typical 
to rotate the whole problem so that the wedge conforms to the bottom boundary. The 
problem is run on the domain [0,3.2] x [0,1]. For the bottom boundary we impose 
the exact solution for x < ^ and we enforce solid wall boundary conditions for x > ^. 
The initial conditions are 




(8.0,8.25,cos(f) ,8.25sin(f) ,116.5) if a; < §, 
(1.4,0,0,1) otherwise. 


In the Cartesian case we run this test using a 960 x 240 cell grid matching the 
second highest resolution in 30 . In the unstructured grid case we use an example 


with 204,479 cells. Both examples are run with the 10-stage SSPRK4 method with a 
CFL constant of 0.08. This smaller CFL number was chosen to guarantee positivity 
of the solution. 

This problem is challenging given that there is a complicated strong-shock struc¬ 
ture, as well as a region where there is vortical flow. It can be seen in Figure]^ that 
both limiters essentially succeed in capturing this behavior. The unstructured grid 
example appears much better than the Cartesian grid example; this is be due to the 
fact that the unstructured grid example makes use of all cells touching a current cell, 
and on the Cartesian grid we only use cells sharing an edge. This was done because 
ideally we do not want to introduce any extra communication in our limiter (by com¬ 
municating across corners). However in the unstructured grid case just limiting using 
cells sharing edges is not good enough and so the unstructured grid limiter had to 
make use of values in cells sharing nodes instead of edges. In Figure]^ we see a close 
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(a) Density for RPl at f = 0.8. One (b) Density for RP3 at f = 0.3. One 
point per cell. point per cell. 




(c) Density for RPl at t = 0.8. Nine (d) Density for RP3 at t = 0.3. Nine 
points per cell. points per cell. 


Fig. 4. 2D Riemann problems. For RPl 30 contours equally spaced between 0.3 and 1.8 are 
plotted. For RP3 30 contours equally spaced between 0.1 and 3.3 are plotted. 


of up the solution in the vortical flow region. Both solutions are certainly forming 
vortices but the unstructured grid again seems to give much better behavior. However 
neither example quite matches the maximum rollup seen in the comparable resolution 
examples in 


30 


8.3. Mach 3 wind tunnel with a step on Cartesian and unstructured 
meshes. We next consider a classical test problem p4,43 that is often overlooked, 


likely due to its difhculty. The problem involves an ideal gas with 7 = 1.4 and a strong 
shock moving at Mach 3 that passes past a forward facing step. The initial conditions 
for this problem are given by = (l.O, 3.0,0.0, 7 “^) . The computational 

domain is the rectangle [0,3] x [0,1] with the region [0.6, 3.0] x [0.0,0.2] cut out. In 
the Cartesian grid case, we use a uniform grid with 480 x 160 elements, and we use 
a total of 90, 342 cells for the unstructured grid case. The Cartesian grid matches 
the moderately refined implementation in [^, however we run with the higher-order 
P 3 elements instead of P 2 elements used in^. Both examples are run with SSPRK4 
using the LLF flux and a CFL of 0.4. 
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Density at t = 0.2 



Fig. 5. Double Mach reflection. For this example we used an a equal to lOOOfe^ ®. We plot a 
total of 30 contours that range equispaced from 1.5 to 22.7. Both plots were done using one point 
per cell. 


This is a complicated problem featuring a long contact discontinuity where, again, 
vortices that should form are observable with high-resolution simulations. An addi¬ 
tional and well-known difficulty is that this example tends to form a spurious entropy 
layer above the step. This results in the creation of a stem-like structure where there 
would normally be a shock reflecting off of the step. This is especially noticeable in 
the Cartesian grid plot. To mitigate this effect, we use a slightly more aggressive im¬ 
plementation of our limiter near top of the step. In place of using point-wise values of 
the primitive variables to estimate local upper and lower bounds, we use cell averages 
of the conserved variables in this edge region only. That is, we replace Mf and mf in 
equations (6.4) and (6.5) by 


Mi := max i qi a{h), max 
jeNr^ 


{<}} 


and rrii := min 


Qi — a(h), min 
jeAfr, 


{<}}, 


where the ^ superscript has been suppressed. Note that this means that in this special 
region we are not using the primitive variables for determining Oi, but actually the 
conserved variables; furthermore, we are using the average values q rather than point 
values at the points defined by Xi ■ Additionally we use a much smaller a value in this 
region to obtain better results. In the majority of the domain we take a = 500/i^'®, 


16 







Density at t = 0.2 


Density at t = 0.2 


0.5 



(a) P 3 on a Cartesian grid 


(b) P 3 unstructured grid 


Fig. 6 . The output from Figure^ zoomed in to show the resolution near the contact discon¬ 
tinuity. This contact discontinuity displays the roll up behavior that is expected in high resolution 
numerical simulations 


but in the small region above the step we take a = The region using the more 

aggressive limiter extends 6 cells above the step in the Cartesian grid case, and until 
y = 0.24 above the step in the unstructured grid case. These values can certainly be 
modified, but we found that if the region is taken to be only one cell width, then the 
entropy layer tends to be much more prominent. Moreover, if the more aggressively 
limited region is too large the shocks in the problem will be smeared out too heavily. 

In Figuresandwe observe that both methods do an excellent job of capturing 
the structure of the contact discontinuity. However it appears that the Cartesian grid 
case suffers from a pronounced spurious entropy layer; although it should be noted 
that the entropy layer is not atypically large (e.g., see [3 51 ). Given the relative 
simplicity of this limiter the results seem very promising. 


9. Conclusions. In this work, we have presented a novel limiter for the discon¬ 
tinuous Galerkin method. It is the first of its kind given its ease of implementation, 
ability to retain high-order accuracy, and robustness for simulation problems with 
strong shocks. The results are achieved without the need for an additional shock- 
detector. The proposed method can be viewed as either an extension of the Barth- 
Jespersen limiter to discontinuous Galerkin methods or as extending the Zhang and 
Shu [46| positivity preserving method to be able to accommodate shocks and not just 
satisfy the positivity principle. Several numerical results were presented, including 
single and multidimensional examples on Cartesian and unstructured grids. Future 
work will involve investigating the ability of this limiter to work on more complicated 
hyperbolic problems such as magnetohydrodynamics (MHD) and three dimensional 
problems with structured and unstructured grids. 
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(a) Cartesian grid 



Fig. 7. Density schlieren plots from the forward facing step problem run on Cartesian and 
unstructured grids. The Cartesian grid used h = and the unstructured grid used 90,342 cells 
(that were uniformly sized). Both were plotted with one point per cell. 
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